clear all
set more off
cd /Users/zimaoxiao/Dropbox/CC_Yield_Predict/replication_package/  // Navigate to replication folder on your own machine

cap mkdir output

qui {
* ============================
* 		   Panel (A)
* ============================
* prepare for mapping
use county_fips CC_tave CC_timeInt30above using data/regression/reg_2020_prediction_3Cstep, clear
bys county_fips: keep if _n==1
merge 1:1 county_fips using data/yield/filter, keep(3)
merge 1:1 county_fips using data/geo/cty_eastern_data, nogenerate
foreach i in CC_tave CC_timeInt30above {
	replace `i' = -9999 if `i' == .
}

* map: 2000-2020 average temperature change (in Celcius)
grmap CC_tave using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -0.15 -0.1 -0.05 0 0.05 0.1 0.3 0.5 0.7 0.9 1.1 1.3) /// 
	fcolor(white ///
	"144 202 249" "187 222 251" "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25") ///
	leg(order(13 "+1.1 -- +1.3" 12 "+0.9 -- +1.1" 11 "+0.7 -- +0.9" ///
	10 "+0.5 -- +0.7" 9 "+0.3 -- +0.5" 8 "+0.1 -- +0.3" 7 "+0.05 -- +0.1" 6 "0 -- +0.05" ///
	5 "-0.05 -- 0" 4 "-0.1 -- -0.05" 3 "-0.15 -- -0.1" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(i) average temperature (­°C)", size(*0.8))
graph save "Graph" "output/FigureS3_A_(i).gph", replace

* map: 2000-2020 >30C temperature bins change
grmap CC_timeInt30above using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -3 -2 -1 0 1 2 3 4 5 6 7 8 9) /// 
	fcolor(white "144 202 249" "187 222 251" "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21") ///
	leg(order(14 "+8 -- +9" 13 "+7 -- +8" 12 "+6 -- +7" 11 "+5 -- +6" 10 "+4 -- +5" 9 "+3 -- +4" ///
	8 "+2 -- +3" 7 "+1 -- +2" 6 "0 -- +1" 5 "-1 -- 0" 4 "-2 -- -1" 3 "-3 -- -2" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(ii) temperature bins (>30°C)", size(*0.8))
graph save "Graph" "output/FigureS3_A_(ii).gph", replace

* prepare for mapping
use county_fips CC_gdd10_29 CC_gdd29 using data/regression/reg_2020_prediction_gdd, clear
bys county_fips: keep if _n==1
merge 1:1 county_fips using data/yield/filter, keep(3)
merge 1:1 county_fips using data/geo/cty_eastern_data, nogenerate
foreach i in CC_gdd10_29 CC_gdd29 {
	replace `i' = -9999 if `i' == .
}

* map: 2000-2020 growing degree days change
grmap CC_gdd10_29 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -20 0 20 40 60 80 100 120 140 160 180 200) /// 
	fcolor(white "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21" "191 54 12") ///
	leg(order(13 "+180 -- +200" 12 "+160 -- +180" 11 "+140 -- +160" 10 "+120 -- +140" ///
	9 "+100 -- +120" 8 "+80 -- +100" 7 "+60 -- +80" 6 "+40 -- +60" 5 "+20 -- +40" 4 "0 -- +20" 3 "-20 -- 0" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(iii) growing degree days (10-29°C)", size(*0.8))
graph save "Graph" "output/FigureS3_A_(iii).gph", replace

* map: 2000-2020 harmful degree days change
grmap CC_gdd29 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -15 -10 -5 0 5 10 15 20 25 30 35 40 45) /// 
	fcolor(white "144 202 249" "187 222 251" "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21") ///
	leg(order(14 "+40 -- +45" 13 "+35 -- +40" 12 "+30 -- +35" 11 "+25 -- +30" 10 "+20 -- +25" 9 "+15 -- +20" ///
	8 "+10 -- +15" 7 "+5 -- +10" 6 "0 -- +5" 5 "-5 -- 0" 4 "-10 -- -5" 3 "-15 -- -10" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(iv) harmful degree days (>29°C)", size(*0.8))
graph save "Graph" "output/FigureS3_A_(iv).gph", replace

* ============================
* 		   Panel (B)
* ============================

* prepare for mapping
use county_fips CC_tave_ssp245 CC_timeInt30above_ssp245 using data/regression/reg_2050_prediction_3Cstep, clear
bys county_fips: keep if _n==1
merge 1:1 county_fips using data/yield/filter, keep(3)
merge 1:1 county_fips using data/geo/cty_eastern_data, nogenerate
foreach i in CC_tave_ssp245 CC_timeInt30above_ssp245 {
	replace `i' = -9999 if `i' == .
}
		
* map: 2020-2050 average temperature change (in Celcius)
grmap CC_tave_ssp245 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -0.8 -0.4 0 0.4 0.8 1.2 1.6 2 2.4 2.8 3.2) /// 
	fcolor(white ///
	"187 222 251" "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25") ///
	leg(order(12 "+2.8 -- +3.2" 11 "+2.4 -- +2.8" 10 "+2.0 -- +2.4" 9 "+1.6 -- +2.0" ///
	8 "+1.2 -- +1.6" 7 "+0.8 -- +1.2" 6 "+0.4 -- +0.8" 5 "0 -- +0.4" 4 "-0.4 -- 0" 3 "-0.8 -- -0.4" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(i) average temperature (­°C)", size(*0.8))
graph save "Graph" "output/FigureS3_B_(i).gph", replace

* map: 2020-2050 >30C temperature bins change
grmap CC_timeInt30above_ssp245 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -3 0 3 6 9 12 15 18 21) /// 
	fcolor(white "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30") ///
	leg(order(10 "+18 -- +21" 9 "+15 -- +18" 8 "+12 -- +15" 7 "+9 -- +12" 6 "+6 -- +9" ///
	5 "+3 -- +6" 4 "0 -- +3" 3 "-3 -- 0" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(ii) temperature bins (>30°C)", size(*0.8))
graph save "Graph" "output/FigureS3_B_(ii).gph", replace

* prepare for mapping
use county_fips CC_gdd10_29_ssp245 CC_gdd29_ssp245 using data/regression/reg_2050_prediction_gdd, clear
bys county_fips: keep if _n==1
merge 1:1 county_fips using data/yield/filter, keep(3)
merge 1:1 county_fips using data/geo/cty_eastern_data, nogenerate
foreach i in CC_gdd10_29_ssp245 CC_gdd29_ssp245 {
	replace `i' = -9999 if `i' == .
}

* map: 2020-2050 growing degree days change
grmap CC_gdd10_29_ssp245 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -100 0 40 80 120 160 200 240 280 320 360 400) /// 
	fcolor(white "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21" "191 54 12") ///
	leg(order(13 "+360 -- +400" 12 "+320 -- +360" 11 "+280 -- +320" 10 "+240 -- +280" ///
	9 "+200 -- +240" 8 "+160 -- +200" 7 "+120 -- +160" 6 "+80 -- +120" 5 "+40 -- +80" 4 "0 -- +40" 3 "-100 -- 0" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(iii) growing degree days (10-29°C)", size(*0.8))
graph save "Graph" "output/FigureS3_B_(iii).gph", replace

* map: 2020-2050 harmful degree days change
grmap CC_gdd29_ssp245 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -20 0 20 40 60 80 100 120 140 160 180 200) /// 
	fcolor(white "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21" "191 54 12") ///
	leg(order(13 "+180 -- +200" 12 "+160 -- +180" 11 "+140 -- +160" 10 "+120 -- +140" ///
	9 "+100 -- +120" 8 "+80 -- +100" 7 "+60 -- +80" 6 "+40 -- +60" 5 "+20 -- +40" 4 "0 -- +20" 3 "-20 -- 0" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(iv) harmful degree days (>29°C)", size(*0.8))
graph save "Graph" "output/FigureS3_B_(iv).gph", replace

* ============================
* 		   Panel (C)
* ============================

* prepare for mapping
use county_fips CC_tave_ssp585 CC_timeInt30above_ssp585 using data/regression/reg_2050_prediction_3Cstep, clear
bys county_fips: keep if _n==1
merge 1:1 county_fips using data/yield/filter, keep(3)
merge 1:1 county_fips using data/geo/cty_eastern_data, nogenerate
foreach i in CC_tave_ssp585 CC_timeInt30above_ssp585 {
	replace `i' = -9999 if `i' == .
}
		
* map: 2020-2050 average temperature change (in Celcius)
grmap CC_tave_ssp585 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -0.8 -0.4 0 0.4 0.8 1.2 1.6 2 2.4 2.8 3.2 3.6) /// 
	fcolor(white ///
	"187 222 251" "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21") ///
	leg(order(13 "+3.2 -- +3.6" 12 "+2.8 -- +3.2" 11 "+2.4 -- +2.8" 10 "+2.0 -- +2.4" 9 "+1.6 -- +2.0" ///
	8 "+1.2 -- +1.6" 7 "+0.8 -- +1.2" 6 "+0.4 -- +0.8" 5 "0 -- +0.4" 4 "-0.4 -- 0" 3 "-0.8 -- -0.4" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(i) average temperature (­°C)", size(*0.8))
graph save "Graph" "output/FigureS3_C_(i).gph", replace

* map: 2020-2050 >30C temperature bins change
grmap CC_timeInt30above_ssp585 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -3 0 3 6 9 12 15 18 21 24) /// 
	fcolor(white "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25") ///
	leg(order(11 "+21 -- +24" 10 "+18 -- +21" 9 "+15 -- +18" 8 "+12 -- +15" 7 "+9 -- +12" 6 "+6 -- +9" ///
	5 "+3 -- +6" 4 "0 -- +3" 3 "-3 -- 0" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(ii) temperature bins (>30°C)", size(*0.8))
graph save "Graph" "output/FigureS3_C_(ii).gph", replace

* prepare for mapping
use county_fips CC_gdd10_29_ssp585 CC_gdd29_ssp585 using data/regression/reg_2050_prediction_gdd, clear
bys county_fips: keep if _n==1
merge 1:1 county_fips using data/yield/filter, keep(3)
merge 1:1 county_fips using data/geo/cty_eastern_data, nogenerate
foreach i in CC_gdd10_29_ssp585 CC_gdd29_ssp585 {
	replace `i' = -9999 if `i' == .
}

* map: 2020-2050 growing degree days change
grmap CC_gdd10_29_ssp585 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 -100 0 50 100 150 200 250 300 350 400 450) /// 
	fcolor(white "227 242 253" ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21") ///
	leg(order(12 "+400 -- +450" 11 "+350 -- +400" 10 "+300 -- +350" 9 "+250 -- +300" 8 "+200 -- +250" ///
	7 "+150 -- +200" 6 "+100 -- +150" 5 "+50 -- +100" 4 "0 -- +50" 3 "-100 -- 0" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(iii) growing degree days (10-29°C)", size(*0.8))
graph save "Graph" "output/FigureS3_C_(iii).gph", replace

* map: 2020-2050 harmful degree days change
grmap CC_gdd29_ssp585 using data/geo/US_County_LowRes_2013coord ///
	, id(_ID) osize(vvthin ...) ocolor("189 189 189" ...) ///
	clmethod(custom) clbreaks(-9999 0 20 40 60 80 100 120 140 160 180 200 220) /// 
	fcolor(white ///
	"251 233 231" "255 204 188" "255 171 145" "255 138 101" "255 112 67" "255 87 34" "244 81 30" "230 74 25" "216 67 21" "191 54 12" "177 0 38") ///
	leg(order(13 "200 -- 220" 12 "180 -- 200" 11 "160 -- 180" 10 "140 -- 160" 9 "120 -- 140" ///
	8 "100 -- 120" 7 "80 -- 100" 6 "60 -- 80" 5 "40 -- 60" 4 "20 -- 40" 3 "0 -- 20" ///
	2 "Not in sample")) ///			
	line(data("data/geo/st_eastern_coord") size(vthin ...)) ///
	graphr(margin(zero)) ///
	legend(size(*1.2) position(3) ring(1)) ///
	title("(iv) harmful degree days (>29°C)", size(*0.8))
graph save "Graph" "output/FigureS3_C_(iv).gph", replace

* ============================
*   Combine to get Figure S3
* ============================
gr combine ///
output/FigureS3_A_(i).gph ///
output/FigureS3_A_(ii).gph ///
output/FigureS3_A_(iii).gph ///
output/FigureS3_A_(iv).gph ///
output/FigureS3_B_(i).gph ///
output/FigureS3_B_(ii).gph ///
output/FigureS3_B_(iii).gph ///
output/FigureS3_B_(iv).gph ///
output/FigureS3_C_(i).gph ///
output/FigureS3_C_(ii).gph ///
output/FigureS3_C_(iii).gph ///
output/FigureS3_C_(iv).gph ///
, row(3) col(4)
noi graph export "output/FigureS3.png", as(png) name("Graph")

erase "output/FigureS3_A_(i).gph"
erase "output/FigureS3_A_(ii).gph"
erase "output/FigureS3_A_(iii).gph"
erase "output/FigureS3_A_(iv).gph"
erase "output/FigureS3_B_(i).gph"
erase "output/FigureS3_B_(ii).gph"
erase "output/FigureS3_B_(iii).gph"
erase "output/FigureS3_B_(iv).gph"
erase "output/FigureS3_C_(i).gph"
erase "output/FigureS3_C_(ii).gph"
erase "output/FigureS3_C_(iii).gph"
erase "output/FigureS3_C_(iv).gph"

}

*** EOF
